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Abstract 

The evolution of Bose-Einstein condensates is usually described by the fa- 
mous time-dependent Gross-Pitaevskii equation, which assumes all bosons to 
reside in a single time-dependent orbital. In the present work we address 
the evolution of fragmented condensates, for which two (or more) orbitals 
are occupied, and derive a corresponding time-dependent multi-orbital mean- 
field theory. We call our theory TDMF(n), where n stands for the number of 
evolving fragments. Working equations for a general two-body interaction be- 
tween the bosons are explicitly presented along with an illustrative numerical 
example. 
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I. INTRODUCTION 



The properties and evolution of Bose-Einstein condensates have been extensively ex- 
plored in the community by employing the famous time-independent and time-dependent 
Gross-Pitaevskii equation, for reviews see [1,2] and for individual applications, e.g., 
Refs. [3-9]. 

Gross-Pitaevskii theory is an excellent mean-field for weakly-interacting bosons when- 
ever a single macroscopic one-particle wavefunction is sufficient to describe the reality. By 
definition, Gross-Pitaevskii theory cannot describe fragmentation of condensates for which 
two (or more) onc-particle wavefunctions are macroscopically occupied. Recently, we have 
developed a time-independent multi-orbital mean-field to describe static properties of con- 
densates [10,11], thus generalizing the (one-orbital) Gross-Pitaevskii mean-field. For other 
approaches to fragmentation of Bose-Einstein condensates, sec, e.g., Refs. [12-14]. Utilizing 
the multi-orbital mean-field has enabled us to find new phenomena associated with frag- 
mentation as well as fermionization of bosonic atoms in traps and optical lattices [15-18]. 
Motivated by this wealth of phenomena in the stationary case, it is the purpose of this 
work to derive a time- dependent multi-orbital mean-field theory, thus generalizing the time- 
dependent Gross-Pitaevskii equation. This will enable one to study dynamical properties of 
fragmented condensates evolving in time. 

II. PRELIMINARIES: TIME-DEPENDENT GROSS-PITAEVSKII EQUATION 

FROM VARIATIONAL PRINCIPLE 

The evolution of interacting structureless bosons is governed by the time-dependent 
Schrodinger equation: 

i7$ = ^— , H{r,,r,,...,rj,) = J2h{ri)+ ^ XoV{r,-rj). (1) 

i=l i>j=l 

Here = 1, is the coordinate of the i-th boson, h{r) is the one-body Hamiltonian contain- 
ing kinetic and potential energy terms, and XoV{ri — r^) describes the pairwise interaction 
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between the i-th and j-th atoms where Aq measures the strength of the interparticle inter- 
action. 

The time-dependent Schrodinger equation (1) cannot be solved analytically, expect for 
a few specific cases only. Thus, approximations are a must. The standard (and simplest) 
mean-field approximation to the time-dependent many-body wavefunction $ assumes all 
bosons to reside in a single time- dependent orbital 4>{r,t), namely, that the many-body 
wavefunction is written as 



This so-called Hartree approximation, as is well known, gives rise to the famous time- 
dependent Gross-Pitaevskii equation [1,2]. For our needs, it would be instructive to show 
how the time-dependent Gross-Pitaevskii equation derives from a variational principle. To 
this end we employ the usual functional action, see, e.g., Refs. [19,20], which reads: 



where the time-dependent Lagrange multiplier ii{t) is introduced to ensure that 4>{r,t) re- 
mains normalized throughout the propagation. To evaluate this action we use the ansatz 
(2) and find: 



$(ri, r2, . . . , rjv, t) = 0(ri, t)(f){r2, t) ■ ■ ■ 0(rjv, t). 



(2) 




(3) 




+ A, 



'0 



N{N-1) 
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j 1 0*(r,t)0*(r',t)y(r-r')0(r,t)0(r',t)(irdr', 



(4) 



By requiring stationarity of the action with respect to variation of 0*(r, t), namely, 



ss 



0, we readily obtain the equation of motion: 



N h{r) + Xo{N -l)J{r,t) 0(r, t) = iiV0(r, t) + /x(t)(/.(r, t). 



(5) 



where the direct, time-dependent local potential J{r,t) reads 




(6) 
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and the shorthand notation = ^ has been introduced. 

Using the constraint that (f){r,t) is normahzed we can ehminate the Lagrange multipher 
lj,{t) from Eq. (5). By taking the scalar product of 4>*{r,t) with (5), the resulting jj,{t) takes 
on the form, 



f^{t)^Nl <j>*{r,t) 



h{r)-i^^+Xo{N-l)J{r,t) 



0(r, t)dr. 



(7) 



Substituting Eq. (7) into Eq. (5) and noticing the identity 



h{r)-f^^+Xo{N-l)J{r,t) 



0(r,t)-#</,(r,t) = 



N 



d 



h{r)-i- + Xo{N-l)J{r,t) 



(8) 



we arrive immediately at the time-dependent mean-field equation: 

Pz0(r, = P [h{r) + Xo{N - 1) J(r, t)] 0(r, t), 
P = l-0(r,t) j (p*{r,t)dr. 



(9) 



Examining Eq. (9) we see that eliminating the Lagrange multiplier fi(t) has emerged as a 
projection operator P onto the subspace spanned by 0(r, t). This projection appears both on 
the left- and right-hand-sides of (9), making it a somewhat cumbersome integro-differential 
non-linear equation. To simplify its appearance we choose the relation 



0*(r,t)0(r,t)dr = 



(10) 



at any time. This condition is equivalent to making the assignment of the time-dependent 
phase 0(r,t) exp|-|- / (^(f)\<j)^ dt^ (f){r,t) in Eq. (9). With this condition, the influence 
of the projection on the left-hand-side simphfies, P0(r, = 0(r, i), and (9) takes on the 
appealing from: 



t<P{r,t) = P 
P = 1 - 



A(r) + Ao(iV-l)J(r,t)] 0(r,t), 
cf>{r,t) J f{r,t)dr. 



(11) 
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The P remaining on the right-hand-side of Eq. (11) makes it clear that the condition (10) is 
indeed satisfied at any time throughout the propagation of (/)(r, t). In practice, the meaning 
of this condition is that the temporal change of 0(r, t) is always orthogonal to 0(r, t) itself, 
which can be exploited to maintain accurate propagation results at lower computational 
costs. Additionally with condition (10), fi{t) in Eq. (7) takes now an appealing form which 
can be interpreted as {N times) the time-dependent chemical potential of the condensate. 

To complete this section we would like to show that the three equations, Eq. (9), Eq. (11) 
and the equation obtained by 'omitting' the projector P completely are fully equivalent to 
each other. To this end, we notice that the Lagrange multiplier /^(t) can be 'absorbed' 
into 0(r, t) by making the assignment 0(r, t) — > exp {i/N J ij,{t)dt} 0(r, t). This immediately 
results in the time-dependent Hartree mean-field equation 



In other words, in the case of a single orbital (p, one can either use the form (9), the form 
(11) or the form (12). 

Finally, to arrive at the time-dependent Gross-Pitaevskii equation in its usual appearance 
we take XoV{r — r') = Xo5{r — r') and set Aq proportional to the s-wave scattering length in 
any of the above equivalent forms (9), (11) or (12). Consequently, the direct potential (6) 
simplifies, J{r,t) — \(f){r,t)f, leading to the familiar result [1,2]: 



In what follows we generalize the time-dependent Gross-Pitaevskii equation, which is 
only valid for non-fragmented condensates, by constructing a multi-orbital time-dependent 
mean-field for fragmented condensates. We will generalize the forms (9) and (11) of the time- 
dependent Gross-Pitaevskii equation where in the latter case constraints such as (10) would 
make the interpretation of the results as well as the numerical implementation particularly 
useful. 



i^{r, t) = \h{r) + Xo{N - 1) J(r, t)] cf>{r, t). 



(12) 



z^{r,t) = \h{r) + Xo{N-l) \(P{r,t)f] 0(r,i). 



(13) 
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III. TIME-DEPENDENT MULTI- ORBITAL MEAN-FIELD EQUATIONS FOR 

BOSONS 



Generalizing the stationary multi-orbital mean-field ansatz [10,11], the most general 
time-dependent mean-field wavefunction for N interacting bosons is the single configuration 
time-dependent wavefunction 



$(ri, r2, . . . , r^v, t) = 50i(ri, t)02(r2, t) • ■ ■ 0iv(rjv, t), 



(14) 



where S is the symmetrization operator. Of course, for bosons not all 0fc(r, t) have to be 
different functions. For instance, in the Hartree (or Gross-Pitaevsku) ansatz (2) all (j)k{r,t) 
are taken to be equal to one another. Generally, however, we may take ni bosons to reside 
in 01 (r, t), 712 bosons to reside in a different orbital 02 (r, t), and so on, as long as the set of 
occupations {n^} satisfies ni+n2 + - ■ ■ + norb = N and, of course, the number Uorb oi different 
orbitals satisfies 1 < norb < N. What is important to note is that the different orbitals 
0fe(r, i) within ansatz (14) are propagated variationally in time and remain normalized and 
orthogonal to one another at any time. 

To proceed for any set of Uorb orbitals and a corresponding set of occupations {nk}, we 
substitute the ansatz (14) for $ into the usual functional action which now reads: 



Tlovb 



/ k,j 



(15) 



where the time-dependent Lagrange multipliers i^kjit) are introduced to ensure that the 
time-dependent orbitals ^^(r, t) remain orthonormal throughout the propagation. The ex- 
pectation value in the action takes on the appearance 



d 

Tiorb 

k 



1. (■^\ ^X!^T/ ^1V^\ T/ 
i>'kk ^ \ '^] + '^O 7) ^kkkk + 9 /^0niVkl[kl] 

\ ^'^ / kk ^ ^ If^k 



(16) 



where the one- and two-body matrix elements are given by 
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Vmi = jj 4>l{r)(t>]{v')V{v - r')4>,{v)(t>i{r')drdv\ 



(17) 



Note the plus sign in Vhj\^qi] which is due to the symmetrization operator. 

By requiring stationarity of the action with respect to variation of all 0^(r, t), namely, 
^ = 0, A; = 1, . . . , Uorh, we obtain a set of Uorh equations of motion for the multi-orbital 

mean-field ansatz (14): 



d_ 

di 



h-i^ + Ao(nfc -l)Jk+^ XqUi {Ji + Kj^ 

l^k 



Wk) 



X] /^fej W \^3) ' A; = 1, . . . , Uorb, 



(18) 



where the direct (local) and exchange (non-local) time-dependent potentials are given by 



Ji{v,t)= j ct>:{v',t)V{v-v')ct>i{v\t)dv', 

Ki{r,t) = J <P1{r',t)V{r-r')rrr4i{r',t)dr', (19) 

and Vrr' permutes the r and r' coordinates of two bosons appearing to the right of it. 

Using the constraints that the 0fe(r, i) are orthonormal functions we can eliminate the 
Lagrange multipliers f^kj{t) from Eq. (18). By taking the scalar product of {(f)j\ with (18), 
the resulting ^Xkjif) take on the form, 

d 



\ / jk Ij^k 



(20) 



Substituting Eq. (20) into Eq. (18) and noticing the identities 



h-i-K; + Mn-k -l)Jk+Yl Ui + Ki) 
^ i^k 

d 
di 



rik J 



i-EI'/'.)(<^.l 



Tlorb 



h-i— + Xo{nk - l)Jk + J2 ^oni (ji + Ki) 

i^k 



(21) 



for all k, we arrive immediately at the time-dependent multi-orbital mean-field equations, 

k — 1 , . . . , Tlorb • 



Pi 



h + Ao(nfe -l)Jk+^ XqUi {Ji + Ki^ 

l^k 



P = l-El<^i)(<^.-|- (22) 

Examining Eq. (22) we see that eliminating the Lagrange multipliers i^kjit) has emerged 
as a projection operator P onto the subspace spanned by the 0fc(r, t). This projection 
appears both on the left- and right-hand-sides of (22), making it a coupled system of integro- 
differential non-linear equations. 

The equation set (22) is the time-dependent multi-orbital mean-field which general- 
izes the one-orbital time-dependent mean-field equation (9). We call our theory in short 
TDMF(n), where n — riarb stands for the number of evolving fragments. Similarly to the 
stationary case [10,11], the multi-orbital time-dependent theory boils down to the standard 
time-dependent mean-field for Uorh = 1- In contrast to the stationary case, where the occu- 
pations rik of the fragments are also determined variationally, the are determined by the 
initial condition, say at t = 0, and remain unaltered during the propagation. To enable the 
variation of the occupations in time, a much more involved many-body theory is needed. 
Such a theory is possible [21]. 

The multi-orbital coupled system (22) is more involved than its one-orbital predecessor 
Eq. (9). Naturally and similarly to the previous section, we would like to make its appearance 
simpler. We can choose the relations 

{(t>k\(t>k) =0, /C = 1, . . . , Uorb (23) 

at any time, which is equivalent to making the assignments of the time-dependent phases 
(f)k{r,t) exp|-|- / (^(l)k\(j)k) dt^ (j)k{r,t) in Eq. (22). With these conditions, the infiuence 
of the projectors on the left-hand-side slightly simplifies, P = Pk where Pk = 
1 — J2]^k Additionally, the diagonal Lagrange multiplier fj.kk{t), see Eq. (20), can 

now be interpreted as (n^ times) the time-dependent chemical potential of the k-th fragment. 



thus generahzing the time- independent quantities [10,16]. However, since the off-diagonal 
terms (^(t)k\(t>j) can a priori not be absorbed in a similar manner by time-dependent 

phases, the system of coupled TDMF(n) equations, combining Eqs. (22) and (23), remains 
involved. 

To proceed let us examine a few properties of the time-dependent multi-orbital mean- 
field for bosons, TDMF(n). First, it is clear that an initially normalized permanent $ 
remains normalized throughout the propagation since the orbit als 0fe(r, i) remain orthonor- 
mal functions at any time. Second, the expectation value of the Hamiltonian with respect 
to $ should be constant if the many-body Hamiltonian H is time-independent. To examine 
this, we calculate the time derivative of Taking the scalar product of (j)k{r,t) 

with Eq. (18), summing up over k and adding the result to its complex conjugate we readily 
obtain: 



H 



"^orb 



hk + ^fcfe + Ao(nfc - 1) {Vj^j^,^^^ + V^^j^^^j^) + 

(24) 



li^k 

In (24) we use the shorthand notation for the index k to indicate expectation values with 
respect to 0fe(r, t), e.g., hj^f, — (^(j)k\h\4>k) , etc. From Eq. (24) we see that the energy would be 
conserved if the matrix of Lagrange multiplier is Hermitian, f^kj{t) = l^'jki't), simply because 

{4>k\4>j) + {(pk\<Pj) =0, k,j = l,...,norb (25) 

hold due to the orthonormality constraints. Whether throughout the time propagation of 
Eq. (22) the matrix of Lagrange multipliers iJ>kj{t) can be kept Hermitian at all times is a 
matter of further studies. There is an alternative to this Hermiticity condition, in the form 
of solutions to Eq. (22) which obey the constraints 



0fe|0i)=O, k,j = l,...,norb (26) 



at any time. When these additional constraints are explicitly invoked, the energy is readily 



The P remaining on the right-hand-side of Eq. (27) makes it clear that the constraints 
(26) are indeed satisfied at any time throughout the propagation of (f)k{r,t). In practice, 
the meaning of these constraints is that the temporal changes of the 0fe(r, i) are always 
orthogonal to the ^^(r, i) themselves, a property which makes the time propagation of 
Eq. (27) robust and stable and which can thus be exploited to maintain accurate propagation 
results at lower computational costs. 

To complete our derivation we would like to pose the question whether, similarly to the 
property of the time-dependent Hartree mean-field discussed in the previous section, see 
Eq. (12), the projector P may totally be 'omitted' in the TDMF(n) theory. Examining 
Eq. (18) we note that the diagonal Lagrange multiphers Hkkit) can indeed be 'absorbed' 
into the orbitals 0fc(i", t) by making the assignments (/>A;(r, t) —>■ exp {i/uk J fikkit)dt} (pkir, t). 
However, in contrast to this situation the off diagonal Lagrange multipliers /ifej(i), k ^ j, 
which ensures orthogonality of the orbitals at any time, cannot in general be 'absorbed' 
into the evolution of the 0fc(r,t), namely that the form (12) does not generally exist in 
the TDMF(n) theory. Therefore and in view of all the above, it is the version (27) of the 
TDMF(n) which will be implemented and employed below in our numerical example. 

Finally, we mention that the presently developed time-dependent multi-orbital theory 
and its stationary version are intimately connected via the so-called imaginary-time propa- 
gation technique. It is well known that the ground state of the stationary Gross-Pitaevskii 
equation can be obtained by propagating in imaginary time the time-dependent Gross- 
Pitaevskii equation, see, e.g., [7]. In our case, by setting t —it, the left-hand-sides of (22) 



conserved, see Eq. (24). The influence of the projection on ^k/ simph 
P = 4>k): and thus (22) takes on the appealing from. A; = 1, . . . , Horb- 




considerably. 




(27) 
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and (27) decay to zero in time. Then, by using the identity (21) and the expression (20) for 
the Lagrange multiphers, the right-hand-sides of (22) and (27) boil down to the stationary 
multi-orbital mean-field theory [10,11]. This establishes a natural connection between the 
time-dependent and time-independent multi-orbital mean-fields for bosons. 



Having derived the time-dependent multi-orbital mean-field TDMF(n) we wish to apply 
it to a specific problem of evolution of fragmented condensates. Of course, we have to chose 
a specific shape for the interparticle interaction and we do so by taking the popular contact 
interaction, Ao^(rj — r^) = Ao5(rj — r^), see Refs. [1,2]. The resulting set of Horb equations 
(27) simplifies considerably and is given by. A; = 1, ... , Uorb' 



The factor of 2 in the last term of the square braces comes from the fact that both the 
exchange and direct potentials are now local potentials which are equal to one another, 
namely, Ji{r,t) — Ki{r,t) — |0;(r, i)|^. Eq. (28) is the generalization of the time-dependent 
Gross-Pitaevskii equation (13) and will be applied below for a one-dimensional problem of 
a two-fold fragmented condensate evolving in time, i.e., for Ugrb = 2. 

As an illustrative numerical example we consider = 1000 interacting bosons in a 
double-well potential. The numerical implementation of TDMF(2) employs the discrete vari- 
able representation (DVR) method [23] and Adams-Bashforth-Moulton predictor-corrector 
integrator [24]. It is convenient to work in dimensionless units in which the one-body Hamil- 
tonian reads h{x) — —\-^ + y{x)- As the trapping potential we chose the symmetric po- 
tential V{x) = O.OSa;^ -|- 10/v^27re~^^/^. The interparticle interaction strength is Aq = 0.1. 
Initially, the system is prepared in the symmetric double- well V{x) in a two- fold fragmented 



IV. ILLUSTRATIVE NUMERICAL EXAMPLE 



i|0fe) = P /i + Ao(nfe-l)|0ifc|'+$^2Aon,|0,p , 




(28) 
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(or, Fock) state $ = \ni,n2), with ni—n2 — N/2 bosons per fragment. The two orthonor- 
mal orbitals (/>i(a;, t = 0) and 02 (a;, t = 0) in which the fragments initially reside are localized 
in the left and right wells of V{x). At time t = 0, the symmetric potential trap V{x) is 
suddenly translated to V{x + 2) and the system is subsequently allowed to propagate in 
time in the new shifted potential V{x + 2) according to the TDMF equation (28). Since the 
system is not in a stationary state any more, the density 



p{x,t) = J2^k\Mx:t)f = - \\Mx,t)f + \Mx,t)f (29) 

k 



evolves in time. 

In Fig. 1 we record four snapshots of the time-dependent density p{x, t). As can be seen 
in the figure, the density profile as a whole oscillates from side to side, which is expected 
from this scenario where the trap has been shifted at i = 0. Atop this behavior of p(x, i), we 
observe that the initial sharp dip of the fragmented-state density at i = nearly vanishes at 
t = 0.9, then it is restored at t = 5.5, and again it nearly vanishes at t = 6.7. Furthermore, 
a few density wiggles which vary in time accommodate p(x, t) in addition to the above 
features, all demonstrating that the time-evolution of fragmented condensates is both an 
intricate and appealing dynamical problem which awaits and invites ample investigations. 



V. SUMMARY 

The evolution of Bose- Einstein condensates has been extensively studied by the well- 
known time-dependent Gross-Pitaevskii equation, which assumes all bosons to reside in 
the same time-dependent orbital. In this work we addressed the evolution of fragmented 
condensates, for which two (or more) orbitals are occupied, and derived a corresponding 
time-dependent multi-orbital mean-field theory, thus generalizing the (one-orbital) time- 
dependent Gross-Pitaevskii mean-field. We call our theory TDMF(n), where n — Uorb 
stands for the number of evolving fragments (orbitals). 

In the TDMF, the ansatz for the many-body wavefunction is taken as a permanent 
made of Uorb time-dependent orthonormal orbitals, whose occupations are determined, say, 
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by the initial conditions. The evolution of the many-body wavefunction is then determined 
variationally by utilizing the usual functional action. This results in a set of Ugrb coupled 
time-dependent equations for the evolution of the Ugrb orbitals (fragments) in time. Working 
equations for a general two-body interaction between the bosons are explicitly presented. 
By considering their propagation in imaginary time, TDMF is shown to naturally relate to 
the recently developed and successfully employed stationary multi-orbital mean-field theory. 

Finally, an illustrative numerical example of the evolution of a two-fold fragmented con- 
densate in a one-dimensional double-well trap is provided, demonstrating an intricate dy- 
namics of the density as time progresses. This is the tip of the iceberg of dynamical properties 
of fragmented Bose-Einstein condensates which await many more investigations. 
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FIGURES 




FIG. 1. (Color online) Time evolution of fragmented condensate made of = 1000 bosons 
in a double-well potential. A two-fold fragmented state $ = \N/2^N/2) is initially prepared in a 
symmetric double- well trap V{x). At time t = 0, the trap is suddenly translated a bit to the left, 
namely, V{x) V{x + 2), and the fragmented condensate is allowed to propagate in time according 
to the TDMF equation (28). Shown are four time snapshots of the density (blue curves), Eq. (29). 
It exhibits an intricate record of overall oscillations from side to side, evolution of the initial, central 
dip, and appearance of additional density wiggles which varies in time. For convenience, the trap 
potential V{x + 2) (black curves) is shown, scaled by 2 and shifted vertically such that V = at 
its minima. 



16 



